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Abstract 

We consider a problem of identification of point sources in time de- 
pendent advection-diffusion systems with a non-linear reaction term. The 
linear counterpart of the problem in question can be reduced to solving 
a system of non-linear algebraic equations via the use of adjoint equa- 
tions. We extend this approach by constructing an algorithm that solves 
the problem iteratively to account for the non-linearity of the reaction 
term. We study the question of improving the quality of source identifica- 
tion by adding more measurements adaptively using the solution obtained 
previously with a smaller number of measurements. 

1 Introduction 

We are interested in a problem of identification of point sources in non-linear 
time dependent advection-reaction-diffusion systems from a sparse set of mea- 
surements. Here by sparse we mean a small number of spatially separated 
measurements. This work is motivated by applications in atmospheric studies 
where one would like to localize a release of an airborne contaminant. A pos- 
sible model for such problem is a linear parabolic system with a known first 
order advection term and point sources [1 . However, for more realistic mod- 
eling of the processes in the atmosphere one need to consider a system with 
multiple chemical species that react with each other. In some cases it may even 
be beneficial to make measurements not of the concentration of the contami- 
nant itself, but of the products of its reactions with the other species in the 
atmosphere. This leads to studying not just a single parabolic equation, but 
a system of such equations. Moreover, an accurate modeling of the chemical 
reactions between the different chemical species requires the use of non-linear 
reaction terms [TOl [TTJ [15] of large magnitudes that lead to very stiff systems. 
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To our knowledge this is the first study of the source identification problem for 
non-linear systems. 

To solve the source identification problem with sparse measurements one 
needs to assume some sparsity of the unknown source term as well. Here we 
assume that both sources and measurements are point-like. Under sparsity 
assumptions in the linear case the source identification problem may be reduced 
to solving a system of algebraic equations obtained by employing the relation 
between the forward model and its adjoint ^i8j. The adjoint problem solution 
is not coupled to the (unknown) forward problem solution, which makes the 
problem much easier to solve numerically compared to general PDE constrained 
optimization problems that arise if no sparsity constraints are used. What 
complicates the non-linear case is that the forward and adjoint solutions are no 
longer uncoupled. In this work we propose a computationally efficient iterative 
procedure that resolves this coupling and solves the source identification problem 
simultaneously. Note that one can try to exploit the sparse nature of sources and 
measurements to use the ideas from compressed sensing [H |5l [16] to recover the 
sources p!2|. This approach can be beneficial if the number of point sources in 
the system is large. However, the main idea of compressed sensing of replacing 
I/O optimization with Li optimization requires some properties of the forward 
operator (like the restricted isometry property), which the forward parabolic 
operator may not satisfy. Thus, we use a different approach here. 

Another aspect of source identification that we consider here is an efficient 
placement of measurements. In the presence of noise in the measured data 
or some uncertainty in the system's parameters an efficient placement of the 
measurements in the domain of interest may play a crucial role in stable source 
identification. Here we consider both a priori placement of initial measurements, 
when one has no prior knowledge about the possible source distribution, and a 
posteriori placement of additional measurements, when one utilizes the source 
estimate obtained with a fewer measurements to place new ones. The problem 
of efficient positioning of measurements is known in the literature under the 
name experimental design or optimal design of experiments |T3l [9]. Here we 
propose a heuristic for adaptive placement of new measurements based on the 
study of a single source case. It is not optimal in the sense that it relies on 
making redundant measurements, however it is computationally inexpensive 
and it performs well in the numerical experiments that we consider. 

2 Non-linear advection-reaction-diffusion system 
with point sources 

A general parabolic system of equations with n components 

studied here has the form 

Ut = DAu-wVu^Lu^Q{u)u^f, xeVt, tG[0,T], (2.1) 
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for some domain ^1 C and terminal time T > 0. Hereafter bold lowercase 
letters denote vectors and vector-functions and bold uppercase letters denote 
matrices and matrix- functions. Dirichlet and Neumann conditions are specified 
on the corresponding parts of the boundary 



du 

dv 



(2.2) 



and the initial condition is 



n(cc,0) = 0. (2.3) 
The diffusion and advection terms are given in terms of the diagonal matrices 



(2.4) 



is the vec- 
tor advection field. The dot product ii;- V in (2.1) is understood componentwise, 
i.e. 

w • Vn = diag {w • Viii, . . . , ii; • Vun) . (2.5) 

Note that the diffusion and advection terms are linear operators. The only 
source of non-linearity in the system is the reaction term 
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R{u) = Lu + Q{u)u^ 

which we split into the linear L and non-linear Q{u) parts. 
We consider the source terms of the form 



fk{x,t) 



ik+i 



ajhj{t)S{x — 2/-^), /c = 1, 



(2.6) 



(2.7) 



where the time dependent part hj (t) of the source term is either a point source 
6{t — Tj) or an indicator function of some time interval. In the simplest case it is 
an indicator function of [0,T]. The source intensities aj > are assumed to be 
constant in time. The spatial location of j^^ source is G ft. The parameters 
= li < I2 < . . . < Ifi ^ ^n+i = determine the number of sources in each 
component, which is lk-\-i — h- The total number of sources in the system is 
denoted by Ns. 

Existence and uniqueness of solutions of non-linear elliptic and parabolic 
systems can be proved using a fixed point iteration technique [14 . Consider for 
simplicity a scalar elliptic operator A and the equation 



Au + R{u) + f{x) =0, X eft, 



where a non-linear reaction term satisfies the condition 

OR 

du 



(2.8) 



(2.9) 
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Then the iteration 



{A - K)U' 



g+1 



(2.10) 



has a unique fixed point that is a solution of (2.8) [14 . This result can be 
generalized to parabolic non-linear systems ( |2.1[ )-( |2.3l ), however the proof tech- 
nique in [14 is not directly applicable to our case since it relies on sufficient 
regularity of the solutions to elliptic (parabolic) equations, which does not hold 
in the presence of point sources. To avoid theoretical complications we assume 
hereafter that the system ( |2.1| )-( [23| ) has a unique solution that can be obtained 
as a limit u{x^t) = lim u^{x,t) of a fixed point iteration 

q^oo 



{DA 



(2.11) 



where for each q we solve the linear system (2.11 ) with boundary and initial con- 



ditions (2.2)-(2.3) for u^~^^ while keeping the previous iterate fixed, starting 



from n^(ic, t) = 0. 



2.1 Adjoint system and source identification problem 

A straightforward way to formulate the source identification problem is to state 
it as an optimization problem with PDE constraints. However, making ad- 



ditional assumptions on the source term like those in (2.7) makes it possible 



to reduce the source identification problem to solving the system of non-linear 
algebraic equations. These equations arise from the adjoint problem. 
Let us define the inner product for vector-functions u and v by 



u{x, t) • v{x, t)dxdt, 



(2.12) 



where u • v 



is a regular inner product in 



For functions that 



are defined on the boundary we replace Q in (2.12) by dfl^ and when time 
integration is not needed we omit T. 



To define equations adjoint to the non-linear system (2.1) we observe that 



if the value of the term Q{u) is known and fixed for the true solution n, then 
(2.1) is a linear system for u. The system of equations adjoint to that linear 



system is given by 



Vt = DAv + w • ■ 



(2.13) 



It runs backwards in time from t = T to t = and thus a terminal condition 
for vix^T) has to be specified. Note that because the time runs backwards in 
(2.13), the system is well-posed, unlike the backward parabolic system that also 



has a minus sign on the left, but runs forward in time. 



Taking the inner product of (2.1) with v and of (2.13) with u we can apply 



the divergence theorem to obtain the adjoint relation 



(2.14) 
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where the correction term is given by 



c{u,v) = -{u,v)^\\Z,^{v,D--) -(u,D--) (2.15) 



The normal derivative ^ in (2.15) is understood component- wise. Typicahy 
one imposes the boundary and initial conditions on the adjoint solution v to 
make as many terms of c(u, v) zero as possible. In particular, to take care of 



the t = T part of the first term in (2.15) we can set the terminal condition to 



v\^^j. = 0. The second and third terms are usually dealt with by enforcing v 
to be zero on the portion of the boundary where 7^ and vice versa. The 
fourth term typically is zero due to the assumption of divergence free advection 
field w{x). Note that if the advection field is divergence free than the correction 
term only depends on the boundary and initial conditions for u and v that are 
known a priori. 



The source term g in the adjoint system (2.13) is chosen according to the 



measurement setup. Since the source / in (2.7) is determined by many param- 
eters a j , (and possibly also ) , j = 1 , . . . , A/'^ , multiple measurements of u 
are needed in order to identify the source term. We denote by gf*^*^ a source term 
corresponding to the i*^ measurement and by v^*^ the corresponding solution of 



(2.13) with g — g'^'^\ i = 1, . . . , Nm where Nm is the number of measurements. 
A single measurement consists of measuring one component ix^. at location 
either at a time instant Oi or integrating over some time interval (usually the 
whole observation period [0,T]). This leads to of the form 

gf {x, t) = 6j^mAt - - z'), j = 1, . . . , n, z = 1, . . . , 7V^, (2.16) 

for the instantaneous measurement, and 

gf{x,t) = 5j,^Mx-z'), j = l,...,n, i = l,...,N^, (2.17) 

for the measurement integrated in time. If we denote the measured data vector 
by 

d, = (g^'\u) , z = l,...,iV^, (2.18) 
then assuming for simplicity that c(n, v) = using the expression for the source 



(2.7) we can rewrite the adjoint relation (2.14) as 

n ^fc+i 

hj{t)vf{y^,t)dt = d,, z = l,...,7V^. (2.19) 

/C=li = ^fc + l 

Note that the above system of equations is linear in source intensities aj and non- 
linear in the source spatial locations (and also possibly temporal locations 
Tj). In what follows it is convenient to express this fact in matrix- vector form 
as 

V{s)a = d. (2.20) 
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Here we stack all the source intensities in the vector a and all source location 
parameters (including the time location parameters tj) in vector s with = 



1 , . . . , and Ng . If the time dependent part of the source term 



is a known indicator function, then we simply have 

Definition 1 (Source identification problem) 

d corresponding to measurements g^^\ z = 1, . . . , Nrr 



Given the measured data 
find the source intensities 



and source location parameters , j = l^...,Ns, that satisfy the adjoint 



relation (2.19). 



The above definition implies that if the adjoint solutions v^*^ are known, the 
source identification problem is equivalent to solving the system of non-linear 



algebraic equations (2.19). This is indeed the case for the linear system, i.e. 
if Q(u) = 0. However, in the non-linear case the adjoint solutions v^*^ 



are 



implicitly dependent on the forward solution u via the (u) term in (2.13). 
The forward solution in turn depends on the source term /, so there is an implicit 
dependency of the adjoint solution on the source, which must be resolved in 



order to solve (2.19). This is done using an iterative procedure that we present 
next. 



2.2 Forward-adjoint iteration for source identification 

To obtain the source parameters a and s we need to solve the system of al- 



gebraic equations (2.19), which requires the knowledge of the adjoint solutions 
Adjoint solutions satisfy a linear system (2.13) which includes the term 



so we must solve the forward system (2.1) for u with an unknown source 
/. We propose the following iterative procedure to solve the source identifica- 
tion problem that iterates over the solutions of both the forward and adjoint 
problems simultaneously. 

Algorithm 1 (Forward- adjoint iteration) 

1. Obtain an initial guess for the forward solution by solving a linear 
system 

= {DA -W'V^ L)u^ 

with boundary conditions {2.^ and initial conditions \2.^ . 
For g = 1, 2, . . . do 

2. Solve the linear systems for the current estimate of the adjoint solutions 



{DA + ii; • V + + Q^{u'i-^))v 



(2.21) 



for i 



, Nm with the appropriate terminal and boundary conditions. 



3. Form the matrix valued function V^{s) for (2.2(f^ from the current esti- 
mates of the adjoint solutions v^'^^'^. 
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4. Obtain the current estimates and of the source parameters by solving 
iteratively the optimization problem 



d\\l 



(2.22) 



minimize \\V^{s)a 

a,s 

and form the current estimate of the source term f^. 
5. Update the estimate for the forward solution by solving a linear system 



ul = {DA - ii; • V + L + Q(n^-^))n^ + 



(2.23) 



with boundary conditions (2.2) and initial conditions \2. 



Convergence of the algorithm can be thought of in terms of both converg- 
ing to the true forward solution u and converging to the true source term 
/. While we are mainly interested in recovering the source term /, convergence 
of one should imply convergence of the other and vice versa. The main idea 
is that (2.23) with an improving source estimate will behave like a fixed point 
iteration (2.11). Convergence analysis appears to be complicated by the fact 
that iteration (2.23) and optimization (2.22) are coupled. Thus, the proof of 
convergence remains to be a topic of further study. 

Since the residual in the objective in (2.22) is linear in source intensities, we 
can eliminate a from the optimization by taking the least squares solution 

a= {y^{s)V{s)^ V^(s)d. (2.24) 

If we substitute the above expression for a into (2.22) the optimization problem 
can be rewritten as 



maximize d^V(s) {y^{s)V{s)^ V^(s)d. 



(2.25) 



Now the optimization obje ctive only depends on source location parameters s. 
The optimization problem (2.25 ) is constrained by G ^1 x [0, T], j = 1, . . . , Ng. 
While it is possible to use a derivative-based approach to solve it, here we use a 
simple derivative- free search procedure that provides good results numerically 
and does not require any extra work to handle the constraints. The algorithm 
below summarizes the search procedure. 

Algorithm 2 (Derivative-free search) 

1. Choose an initial guess for source location parameters s. 

For p = 1, 2, . . . do 

For j = l,...,Ns do 

2. Freeze all the components of s for k ^ j and compute the 



objective 

J{s) = d^V{s) (v^{s)V{s)y' 
for all possible values of G f2 x [0,T]. 



V^{s)d 



(2.26) 
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3. Update the location of the j source 

= argmax J{s^,...,s^-\r,s^^^,...,s^^). (2.27) 

renx[o,T] 

4. If for all j = 1, . . . , A/'g the changes in step 3 compared to iteration 
p — 1 are small then stop. 

Algorithm [2] has an inner-outer iteration structure. At each outer iteration 
indexed by p the algorithm cycles through all source locations s-^ , j = 1, . . . , A/'^ 
and performs an exhaustive search for each of them while keeping the rest fixed. 
While it may seem as a computationally expensive solution, we should note that 
Algorithm [2] is just a single step in Algorithm [l] and in practice it is the cheapest 
step. Most of the computational time in Algorithm [l] is spent computing the 
adjoint solutions in step[2j so the computational cost of step |4] is negligible. 

Different stopping criteria can be used in step [l] of Algorithm |2j In practice 
since the adjoint systems ( |2.21[ ) are solved on a finite grid, one can use "no 
change from iteration p — 1" as a stopping criterion in step[l] Also, the number 
p of outer iterations of Algorithm [2] can be used as a stopping criterion for the 
iteration indexed by q of Algorithm [l] In particular, if Algorithm [2] terminates 
with p = 1 then we can terminate the Algorithm [l] as well. In practice this 
approach will terminate before the adjoint ( 2.21[ ) and forward ( |2.23[ ) solutions 



fully converge, so the estimate of the source strength ( |2.24 ) might be slightly 



inaccurate due to inaccuracy in the adjoint solutions. However, such approach 
gives quite accurate estimate of the source positions s. Moreover, this saves 
considerable amounts of computation, because the expensive step ( |2.21|) is not 
performed as many times as needed for full convergence of (2.21) and ( |2.23 ). 



Since the Algorithm [2] is used inside the iterations of Algorithm [T] one can 
take as an initial guess for s in step 1 of Algorithm [2] the estimate for s from 
iteration g — 1 of Algorithm [T] Then one has only to determine the initial 
guess for s at the beginning of Algorithm [l] While it is possible to use a 
randomly chosen guess or a guess obtained from some prior knowledge of the 
source position, we propose a systematic way of obtaining the initial guess from 
the measured data d only. It is summarized below. 

Algorithm 3 (Initial guess for source locations) 

1. Given the initial guess from step\^ of Algorithm^ with q = I, assemble 
the matrix assuming that there is only one source present. In this case 
only has one column and depends on only. Thus the optimization 



objective J in (2.26) also depends on only. 



2. Compute the estimate of the first source location as 



5^ = argmax J{r). (2.28) 

renx[o,T] 



For k = 2,...,N, do 
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3. Assemble assuming that there are k sources present. Fix the locations 
of previously determined sources j = 1, . . . , /c — 1 so that the optimiza- 
tion objective J only depends on . 

4- Compute the estimate of the k^^ source location as 

= argmax J{s^,. . . ,s^~^,r). (2.29) 

rGOx[0,T] 

Note that in the case of a single source TVg = 1, there is no need for an initial 



guess since (2.28) is the same as (2.27). In this case we can think of J(r) as an 
imaging functional, which quantifies the likelihood of the source being located 
at point r G Q x [0,T]. With noiseless measurements and exact knowledge of 
the adjoint solutions the true location of the source corresponds to the point 
where the imaging functional attains its maximum. 

2.3 Measurement placement and determining of the num- 
ber of sources 

In this section we study the question of choosing the locations at which mea- 
surements are made, which we hereafter refer to as measurement placement. 



Each source in (2.7) is determined by at least d -\- 1 parameters, which are the 
spatial location coordinates and intensities aj , and possibly also the temporal 
locations r^-, j = 1, . . . , A^g. Thus, in the simplest setting of time-independent 
sources we need at least d -\- 1 measurements per source so that the non-linear 



system (2.20) is formally determined. In practice it is beneficial to have an 



overdetermined system (2.20) since having redundant data makes source detec- 
tion less sensitive to noise. Aside from the measurement noise there is also an 
issue of robustness of optimization Algorithms [T] and [2j In the numerical ex- 
periments we observed that having redundant measurements also increases the 
robustness of optimization, i.e. Algorithms [l]-[3] are less likely to get stuck in 
local minima if more measurements are added. 

The problem of choosing the number and positions of measurements has two 
aspects to it: 

1. Initial placement of measurements before any data is available. 

2. Adding more measurements to the existing setup based on the estimate 
of source locations obtained from the data already measured. 

To formulate a strategy of adaptive measurement placement we study how 
the measurement positions affect the estimates of source location in case of a 
single source. Such strategy may not be optimal, because it does not take into 
account the interactions between multiple sources, but it allows us to come up 
with an easy set of rules that can be applied if adding new measurements is 



relatively inexpensive. As we see in the numerical experiments in Section 4.3 
such a strategy indeed proves itself useful in the case of multiple sources. 
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Figure 1: Sum of level set indicator functions X{y) for the case of one source 
and three measurements. Measurement positions are yellow o, source position 
is yellow x, source estimate is black □. 



We first consider the simplest case of identifying a single source with a known 
intensity. Equations ( 2.2Q| ) then become V{y) = d/a, with both d and a 



scalars. Thus, the source is located at the intersection of level sets of V{y) 
corresponding to the value d/a. The level sets are closed curves relative to 
the domain Q. There exists an analogy to the process of triangulation in radar 
detection, where the corresponding curves are circles. The analogy is exact for a 
linear diffusion equation in R^, for which the level set curves are circles too. To 
illustrate this analogy numerically we consider a problem with one source and 
three measurements in two dimensions (the detailed description of the system is 
given in Section [s]). Let us introduce the indicator functions of C-neighbor hoods 
of d/a level sets 

1' if\V,{y)-d/a\<C 
^^^^^"\0, if \Vk{y)-d/a\>C ' i^-^Uj 



3 

for some C > 0. Given the sum X{y) = ^ Xk{y) we can define the set 
Ss = {y W X{y) — 3}, then the position of the source can be estimated as 

y^ = . (2.31) 

J y 

This is shown in Figure [l] with C = 0.125. 

For the case of exact data we can place the measurements anywhere in the 
domain and still be able to recover the location of the source. However, the 
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Figure 2: Imaging functional J{r) from (2.28) for different measurement posi- 



tions. Left: all measurements upwind. Middle: all measurements downwind. 
Right: mixed measurements (2 downwind and 1 upwind). Measurement posi- 
tions are yellow o, source position is yellow x, estimated source location (max- 
imum of J{r)) is black □. 



presence of noise effectively limits the distance from the measurement to the 
source that allows a stable source identification. This is due to the fact that 
the magnitude of the measured solution u decays quickly away from the source, 
and measuring weak signals is more prone to errors than measuring strong 
signals. Thus, if a priori information about the source locations is not available, 
a reasonable strategy is to distribute the measurements more or less uniformly 
around the domain Q. 

The situation is different when some prior knowledge about source locations 
is a available. One example is when we would like to add new measurements 
adaptively based on the results of source identification with a previously chosen 
smaller number of measurements. Here we assume that we can add new mea- 
surements anywhere in the domain and that it is relatively inexpensive to do so. 
This leads us to a strategy of adding a few new measurements for each previ- 
ously identified source. We also assume that the sources and measurements are 
active for all times t G [0, T], so only the spatial placement of the measurements 
is considered. For a source with unknown intensity we add d-\-l measurements, 
where we use d = 2 dimensions for the convenience of visualization. 

In order to distribute the newly added measurements around the previously 
estimated source locations we study how the distribution of measurements af- 
fects the source identification in the presence of advection in case of a single 
source. In Figure [2] we plot the imaging functional for the three different mea- 
surement distribution (the details of the numerical setup are given in sections 
|3]and[4|. The advection direction in Figure [2] is from right to left, so we re- 
fer to the measurements to the right of the source as upwind and to the left 
of the source as downwind. The three possible distributions given are for all 
three measurements upwind, all three measurements downwind and a mixed 
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distribution of one measurement upwind and 2 downwind. 

The plots in Figure [2] are for the noiseless data, so the source position is 
recovered exactly (up to the nearest computational grid point). However, there 
is a drastic difference in the behavior of the imaging functional, which allows 
us to identify an optimal placement of measurements. Obviously, having all 
measurements upwind is the worst scenario. Advection propagates the plume 
away from the measurements and makes source identification difficult. This 
is reflected in the imaging functional having a vast plateau which implies the 
lack of discriminatory power of such functional. Ideally an imaging functional 
should have a single concentrated peak at the source location. By placing all 
three measurements downwind the behavior of the imaging functional is much 
improved. The peak is now located on a narrow ridge, so the localization of 
the solution is much better. Finally, we observe that having one measurement 
upwind can further improve the imaging functional since it allows for exclusion 
of a portion of the domain from possible source locations (the imaging functional 
is small around the upwind source). 

Considering the above observations we propose the following procedure for 
adaptive measurement placement. 

Algorithm 4 (Adaptive measurements placement) 

1. Obtain an estimate of source locations y. 

2. Choose a trust radius pr and a reference simplex T with vertices , 
/c = l,...,(i+l. The orientation of the reference simplex is such that one 
vertex lies upwind and d vertices lie downwind from its center ( the center 
of circumscribed sphere). 

For j = l,...,Ns do 

3. Place the center of the reference simplex at . 
For A: = 1, . . . 1 do 

4. Place a new measurement in the direction of the vertex at a dis- 
tance 

p = min {pT.K.n dist{y^,dQ) ,Ky dist{y^, {y' \ i ^ j})) (2.32) 

away from y^ , where the constants hiQ^tvy G (0,1) determine how 
close the new measurements can be placed to the boundary and the 
rest of the sources respectively. 

For z = 1 , . . . , j — 1 do 

5. Place a new measurement between y^ and y^ . 

The choice of a trust radius pr in step 2 should be determined by the noise 
level, i.e. the distance from the source to the measurements for which the source 
can be identified in a stable manner. In two dimensions the algorithm places 
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three new measurements per identified source in a triangular pattern around 
each source so that one measurement is placed upwind and two downwind. 
Relation (2.32) ensures that the new measurements are not placed too close 
to the boundary or to other sources. This helps to separate the sources in 
case they are clustered together. Adding measurements in step 5 also helps 
separating clustered sources. It is possible to adjust the shape of the reference 
simplex T in step 2 and the positioning of the measurements in step 5 to take 
into account the knowledge of the advection field. However, for simplicity in 
the numerical examples in Section [4] we use an equilateral triangle T in step 2 
and z — {y^ + ?/*)/2 in step 5. 

Note that the algorithm [4] is based on the idea of local refinement, i.e. the 
new measurements are placed near the estimated source positions. Such ap- 
proach is in agreement with the successive sampling strategy developed in [12 , 
which gives good results for source identification in Li setting. 

We conclude this section by considering the problem of determining an un- 
known number of sources, the case when Ng is not known a priori. A procedure 
that appears to be both simple and reliable if to start with and estimated num- 
ber of sources TV* = 1 and run Algorithm [l] repeatedly for increasing numbers 

= 2,3, — Note that the optimization problem (2.25) does not impose any 
constraints on the signs of components of a. Thus, for some value of A/"* Algo- 
rithm [T] will compute a solution with aj = for some j in the noiseless case, 
or in the presence of noise aj < e (this includes negative aj) for some small 
e, which should be chosen based on noise level. Once this happens we deter- 
mine the true number of sources as Ng = N* — 1. The choice of the number 
of measurements in the case of unknown Ng can be done in two ways. If an 
upper bound Ng < N^^^-^ is available one may set Nm to the number of mea- 
surements needed to identify N^^^-^ sources stably. Alternatively, one may add 
the measurements adaptively using Algorithm |4] for each value of . In the 
numerical experiments in Section [44] we use the former approach. 



3 Three component chemical system 

In this section we consider a system that we use in the numerical experiments in 
Section [4j We use a simplified, but somewhat realistic three component n = 3 
chemical system that models the chemical processes occuring in the atmosphere 
based on Chapman's cycle [151 [11]. While a realistic atmospheric model may 
contain dozens of reacting species [10], our simple model still captures some 
of the basic features of atmospheric models like polynomial non-linearity and 
stiffness. 



3.1 Forward problem 

The components of the system are the nitric oxide {NO), nitrogen dioxide {NO2) 
and ozone (O3) denoted by ixi, U2 and us respectively. We assume that nitrogen 
dioxide is released at source locations and the concentrations of nitric oxide are 
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measured. A simplified model of chemical reactions in the system is 



7VO + O3 ^ NO2, (3.1) 
NO2 ^ NO + O^, (3.2) 

with rates ki = 1000, k2 = 2000. Let us introduce a scalar reaction term 



r{u) = kiuius - k2U2, 
then the vector reaction term is given by 



R{u) 



where we take 
L 






k2 
















k2 






-r{u) 
r{u) 
-r{u) 



Q{u) 



Lu + Q{u)u, 



-kius 
kiUs 




(3.3) 



(3.4) 






-kiui 



(3.5) 



Note that while the linear part L is defined uniquely, different definitions of the 
quadratic part Q are possible that lead to the same value of the matrix-vector 
product Q{u)u and thus the same reaction term. 

The realistic values of the diffusion constants are ei = 1, e2 = 63 = 5, which 
lead to a rather stiff system of equations due to large contrast between the 
diffusion constants and the reaction rates /ci, k2- 

While our method works in any number of spatial dimensions, in this nu- 
merical example we use d = 2 dimensions for the simplicity of visualization. 
The system is solved in the unit square with circular obstacles 



No 



[0,l]^|U^r,(c,)), 



(3.6) 



where No is a number of obstacles. In the example below we take No = 2. 
Dirichlet conditions are enforced on the outer boundary 



'^ila[o,i]2 — '^2la[o,i]2 — '^3la[o,i]2 — 
and zero Neumann conditions are enforced on the obstacle boundaries 



9uj 
dv 



0, j = l,...,n. A: = l,...,A^o- 



Constant initial conditions are used 

^ilt=o = ^2|t=o = 0' ^slt^o 



1. 



(3.7) 



(3.8) 



(3.9) 
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We assume that at t = all sources go off and remain active for the period 
of time [0,T]. The source term has the form 







J2 ajS{x-y^) 




(3.10) 



so only the source locations G O and the constant source intensities aj > 
are to be determined. 



3.2 Advection field 

A realistic assumption on the advection term is that there exists a preferred 
advection direction wq that does not depend on time. It is also reasonable 
to assume that the advection vector field satisfies non-penetrating boundary 
conditions on the boundaries of the obstacles 



(3.11) 



Let us introduce the advection potential (j) such that 



(3.12) 



Then the condition that the advection vector field is divergence free implies that 
(j) must be harmonic 

A(/) = in (3.13) 
with zero Neumann boundary conditions on the obstacle boundaries 



dv 



0, j = l,...,7V,, 



(3.14) 



aBr.(cj) 

and Neumann conditions enforcing the preferred direction on the outer bound- 



aries 



= (^0 •z^)la[o,i]^ • 



(3.15) 



Advection field used in the numerical examples below corresponds to a pre- 
ferred advection direction wq = (—50,0), i.e. the "wind" blows from right to 
left. 



3.3 Measurements and the adjoint system 

For the three component chemical system we measure the component ui. In 
the numerical results below we consider the cases of measurements integrated in 
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time (sections |4.3| and 4.4) and of meas urements instant in time (Section |4.5[ ) 
The term g in the adjoint system (2.13) for i^^ measurement takes the form 



'S{x - z') 





or 



'5{t 



i)5{,x-z') 





(3.16) 



for integrated or instant measurements respectively. 

While the initial and boundary conditions for the forward system typically 
come from the physical problem, we have a freedom of choosing the terminal and 



boundary conditions for the adjoint system, so that the adjoint relation (2.14) 



is as simple as possible. In particular, we would like the correction term (2.15) 



to be zero. We enforce zero Dirichlet conditions on the outer boundary and zero 
Neumann conditions on the obstacle's boundaries for all three components of 
the adjoint solution v. We also use zero terminal condition v{x^T) = 0^ x G fl. 
From the expression (3.5) for L and Q{u) we obtain the equation for the 



third component of the adjoint solution 



{esA — w V — kiui) V2,. 



(3.17) 



Combined with the terminal and boundary conditions we immediately see that 



v^{x,t)={), xeVt, tG[0,T]. 



(3.18) 



This is enough to make the correction term ( 2.15| ) zero. Indeed, the only com- 
ponent of u and v that has non-zero initial (terminal) or boundary conditions 
is u^. Since ^3 is identically zero it neutralizes non-zero initial and boundary 
conditions for in the first three terms of the c(n, v). There is no contribution 
from the other components of n, and v to the first three terms, so they are iden- 
tically zero. The fourth term is zero since we use a divergence free advection 
field w. Finally, the fifth term on the outer boundary is taken care of because 



la[o,i]' 



= 0. On the boundaries of the obstacles it is zero since w satisfies the 



non-penetrating conditions ( |3.11[ ) there. 

Once we establish that c(n, v) is zero we can write the components of the 
system of equations (2.20) arising from the adjoint relation 



V^k= I v^\y\t)dt, or V,k = v^\y\Tk), (3.19) 
^0 

di= 1 ui{z\t)dt, or di = ui{z\Oi), (3.20) 
Jo 

for integrated or instant measurements respectively, where i = 1, . . . , k = 
1 , . . . , A/'g and 



a = (ai, . . .aivj 



(3.21) 



is the same for both cases. 
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4 Numerical results 



We implement our method of source identification and provide the results of 
the numerical experiments below. The first two sets of experiments use time 



integrated measurements as described in Section 3.3 In these experiments we 



identify time-independent sources in the cases where the number of sources itself 



is known (Section 4.3) or unknown (Section 4.4). In Section 4.3 we also study 



adaptive positioning of measurements and its influence on source identification. 



Finally, in Section 4.5 we provide the numerical results for identification of a 
time dependent source from measurements instant in time. Results from both 
one and two dimensional settings are presented. 

4.1 Linear parabolic solver 

We solve the linear parabolic systems for the forward and adjoint iterations 
using the following numerical schemes. The spatial part is discretized with fi- 
nite differences on a uniform cartesian grid. The Laplacian in the diffusion 
term is discretized using the standard five-point stencil. The advection term is 
discretized using a central difference scheme. The reason for using the central 
difference scheme for the advection term is that we can use the same discretiza- 
tion for the forward and adjoint problem, for which the direction of advection 
is reversed. Note that such discretization may become inaccurate if the advec- 
tion dominates other terms. While there exist more sophisticated and accurate 
numerical schemes for the solution of ad vect ion-diffusion equations, the focus of 
this work is not the numerical solution of the forward problem. The numerical 
scheme described in this section appears to be sufficiently accurate for source 
identification in a three component system described in Section [3j 

To obtain the solution in time we use an exponential integrator. After dis- 
cretizing in space we need to solve the system of ODEs for the forward and 
adjoint problems of the following form 

^, = Eit)^ + Cit). (4.1) 

The dependency of the matrix E on time is due to the fact that the reaction 
term Q depends on the forward solution that is a function of time. If we denote 
the k^^ time step by tk and the size of the step is hk = t/c+i — ^/c, then the 
approximate solution at time step /c + 1 is given by 



.(/c+l) 



exp [E^'^h,) (^[e^'^) ' C^'^ + C^'^) - {e^'^) ' C^'\ (4.2) 



where ^ |(t/,), E^^^ = E{tk) and C^^^ = C(^/e). While each step of this 
method is more computationally expensive than that of traditional time stepping 
methods, it is a lot more accurate allowing us to use a small number of time 



steps. Note that (4.2) requires evaluation of matrix- vector products with matrix 



exponentials. We evaluate these products using an efficient algorithm 
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In order to avoid committing an inverse crime [6^ we use different grid and 
time steps for the forward problem data simulation and for the solution of the 
source identification problem with Algorithm [l] We simulate the data on a 
finer grid with 80 grid nodes in both x and y directions and 30 time steps. In 
the case of time independent source and integrated measurements we perform 
source identification on a grid with 63 grid nodes in both x and y directions and 
19 time steps. 

Note that even without adding artificially generated noise to the data, using 
different (and relatively coarse) grids for the data simulation and source iden- 
tification is equivalent to having some systematic error in the measurements. 
This poses an issue in the case of time dependent source and instant measure- 
ment, since this case is more sensitive to the noise level in the data. In this 
case we use a finer grid for source identification, namely with 73 nodes in both 
directions. Another modification to the solver required in this case is the use of 
non-uniform time stepping. In order to properly resolve the singularity of the 
sources around times r/c, /c = 1, . . . , A/'^ in the forward problem and around 
j = 1, . . . , Nm in the adjoint problems, we refine the time stepping locally. 



4.2 Noise model 

We provide below the results of the numerical experiments for identifying sources 
from noisy measurements. Single source identification with noiseless measure- 
ments can be found in Figure [2] In this section we use a simple noise model 
with multiplicative normally distributed noise. Such model while being easy to 
implement captures a realistic assumption that the noise level can be viewed as 
constant relative to the strength of the signal. 

If we denote the simulated data vector by d, then the noisy data d* is given 

by 

d* = (I + aAr)d, Ar = diag(Xi,...,X^^), (4.3) 

where a is a scaling term and Xj, j = 1,...,A/'^ are independent normally 
distributed random variables with zero mean and unit standard deviation. All 
results are presented for one particular realization of noise, although for different 
realizations the results remain similar, which implies that the source identifica- 
tion Algoritm [l] is relatively stable. 

In the case of time independent sources and integrated measurements we 
take the scaling factor a = 0.05 corresponding to 5% relative noise level. Note 
that according to (|4.3[) the noise is added to the data after the integration in 



(3.20). Adding the noise to ui before the integration in (3.20) would make it 
easier for Algorithm [l] to determine the source, since integration in (3.20) would 
act as a noise canceling filter. In order to stress test our method we add the 
noise after the integration instead. 

The case of time dependent source and instant measurements is more diffi- 
cult, so we reduce the noise level to a = 0.01 for the numerical experiments in 
Section 14.51 
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Table 1: True and reconstructed source intensities aj and locations for the 
case of two sources (5% noise in the data). 



Case 


ai 


^2 






True 


10.00 


7.00 


(0.70, 0.30) 


(0.30, 0.70) 




13.37 


2.27 


(0.70, 0.32) 


(0.29, 0.87) 




9.62 


6.44 


(0.67, 0.30) 


(0.29, 0.70) 



4.3 Identifying multiple sources with adaptive measure- 
ment placement 

In this section we study the identification of a known number of time indepen- 
dent sources from integrated measurements in a three component system from 
Section |3] We consider two cases Ng = 2 and Ng = 3 in Figures [3] and [4] re- 
spectively. To demonstrate the adaptive measurement placement we begin by 
choosing the smallest number of measurements = SNs that yields a formally 
determined system (2.20) (top row of Figures [3]and[4|. Then we add more mea- 



surements according to Algorithm [4] and run the Algorithm [T] again (bottom 
row of Figures [3] and [4]). For the purposes of visualization for the j^^ source we 
fix all other source locations given by Algorithm [l] and evaluate the functional 



( |2.27[ ) for all possible locations of the j^^ source. Obviously, the maximum of the 
functional corresponds to the source location estimated by Algorithm [l] Doing 
so allows us to visualize the sensitivity of the objective function with respect to 
a particular source location and also how it changes when more measurements 
are added adaptively. 

We observe right away in Figure |3] for the case Ng = 2 that the objective 
function ( |2.22[ ) is not convex (the functional J is not concave). Also in Figure [3] 
for Njn = 9 measurements it is clear that the objective can develop narrow val- 
leys (ridges of J) and become multimodal. This makes the source identification 
problem difficult to solve, and we observe that in the presence of noise the esti- 
mated source location may differ from the true one if too few measurements are 
used (top rows of Figures [3] and [4| . On the other hand, if more measurements 
are added adaptively by Algorithm |4j the behavior of the objective improves, 
as can be seen in the bottom rows of Figures [3] and |4] The imaging functionals 
corresponding to each source become much better localized and less multimodal. 
This leads to an easier optimization problem and much better estimates of the 
source locations and intensities. 

In Figures [3] and |4] we only show the estimated locations of the sources. The 
corresponding source intensities (and the numerical values for the locations) 
are given in Tables [T] and [2] for the cases Ng = 2 and Ng = 3 respectively. 



From the presented data we observe that the least squares estimate (2.24) of 
source intensities is quite sensitive to the estimate of the source locations. When 
few measurements are used and the estimates of the locations are not accurate 
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Figure 3: I magin g functionals J{r^y'^) (left column) and J{y^^r) (right col- 
umn) from J^2.27 ) for all r G 1^ evaluated at the final value of {y^^y'^) given by 
Algorithm [1] Top row: initial run with Nm = 6 measurements. Bottom row: 
subsequent run with adaptively added measurements, Nm = 13. True source 
locations are yellow x, measurement locations are yellow o, estimated source 
position (y^,?/^) - maximum of the imaging functional is black □, initial guess 
from Algorithm [3] is black ^. 



enough the estimated intensities differ significantly from the true values. How- 
ever, when more measurements are added adaptively, the intensity estimates 
improve greatly. Note that this limitation of our method comes from the fact 
that in (2.24) we eliminate the source intensities from the optimization variables. 
If we have some a priory knowledge about the intensities (i.e. the bounds) we 
can retain a as an optimization variable in (2.22) and enforce our a priori knowl- 
edge as a constraint. This comes at a price of enlarging the space of optimization 
variables, which makes the optimization problem harder to solve. 
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Table 2: True and reconstructed source intensities aj and locations for the 
case of three sources (5% noise in the data). 



Case 


ai 


a2 


as 






y' 


True 


10.00 


7.00 


5.00 


(0.70, 0.30) 


(0.20, 0.80) 


(0.50, 0.40) 




13.06 


21.42 


4.22 


(0.70, 0.32) 


(0.22, 0.98) 


(0.43, 0.40) 




9.35 


7.26 


5.31 


(0.70, 0.30) 


(0.19, 0.80) 


(0.48, 0.40) 




Figure 4: Imaging functionals J{r,y'^^y^) (le ft col umn), J{y^^r^y^) (middle 
column) and J{y^^y'^^r) (right column) from (2.2_7) for all r e Q evaluated at 



the final value of (y^^y'^^y^) given by Algorithm [l Top row: initial run with 
Nm = 9 measurements. Bottom row: subsequent run with adaptively added 
measurements, Nm = 21. True source locations are yellow x, measurement 
locations are yellow o, estimated source position {y^^y'^^y^) - maximum of the 
imaging functional is black □, initial guess from Algorithm |3] is black ^. 



4.4 Source identification for an unknown number of sources 

Let us now consider source identification in the case when a true value of Ng is 
unknown. Similarly to the previous section we identify time independent sources 
from integrated measurements. To identify the sources and their number we use 
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Figure 5: Unknown number of sources, the true value is Ng = 3. Imaging 
functionals evaluated at the final value of (y^, . . . , y^^) given by Algorithm |l| 
Top row: = 1 (leftmost) and N* = 2 (middle and right). Middle row: 
A^* = 3. Bottom row: = 4. Number of measurements = 20 for all 
A/"*. True source locations are yellow x, measurement locations are yellow o, 
estimated source positions (y^, . . . , y^^ ) - maxima of the corresponding imaging 
functional are black □, initial guesses from Algorithm |3] are black ^. 



the procedure from Section 2.3 To simplify the exposition in this section we 
do not add the measurements adaptively. Instead we use a predetermined large 
number of measurements for all trial values of A^*. The measurements are 
distributed somewhat uniformly in 1], as shown in Figure [5] 

We set A's = 3 and we perform four trials A"* = 1,2,3,4. The results of 
these trials are given in Figure |5] and Table |3] We observe that as we increase 
the the trial number A'* Algorithm [l] starts to "notice" the sources with smaller 
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Table 3: True and reconstructed source intensities aj and locations for the 



case of unknown number of sources (5% noise in the data). 



Case 


ai 


a2 


as 


a4 


y' 


y' 


y' 


y' 


'Itue 

{Ns = 3) 


10.00 


7.00 


5.00 




(0.70, 0.30) 


(0.20, 0.80) 


(0.50, 0.40) 




n: = i 


12.29 








(0.67, 0.35) 








N* = 2 


12.12 


8.13 






(0.67, 0.33) 


(0.19, 0.79) 






=3 


10.44 


6.83 


4.44 




(0.69, 0.30) 


(0.19, 0.80) 


(0.48, 0.43) 




=4 


11.30 


6.42 


4.12 


-0.17 


(0.69, 0.30) 


(0.20, 0.79) 


(0.48, 0.41) 


(0.62, 0.14) 



intensity. At the first step A/"* = 1 is picks the dominant source with ai = 10. 
At the second step it notices the presence of the source with a2 = 7. Note that 
the locations and intensities of the first two sources are not determined exactly. 



because the objective (2.25) is different from the true one unless = Ng. 
However, the estimates of the locations and intensities of the first two sources 
while not exact are quite accurate, as can be observed in the first row of Figure 
|5] and also in Table [3] The order in which the procedure finds the sources may 
not necessarily be determined by their intensity. In this example the method 
gives preference to stronger sources because the measurements are distributed 
more or less uniformly throughout the domain, so the influence of each source 
is mostly determined by its intensity. 

Finally, when we reach A"* = A'^ = 3 the method identifies all three sources 
quite reliably given the level of noise present. As we go one step further A"* = 4 
Algorithm [l] recovers a spurious source with a negative intensity, which is an 
indicator of an overestimation of the number of sources. Thus, we conclude 
that the true sources were recovered in the previous step and the true number 
of sources is A'^ = 3. Note that while a spurious source appears in the case 
A^* = 4, the method gives a good estimate of the locations and intensities of the 
three sources. We observed such behavior for many realizations of the noise, so 
the results presented in Figure [5] and Table [3] are representative of the general 
performance of the method. 

4.5 Time dependent source identification 

In the numerical examples considered above we used time independent sources 
that are active for all t in [0, T]. In this section we apply Algorithm [l] to identify 
a single time dependent source, which is a point source in both space and time. 



For simplicity of visualization we first consider in Section 4.5.1 a problem in one 



spatial dimension. This allows us to plot the imaging functional for all space 



and time locations. In section 4.5.2 we consider the example in two dimensions 



for the three component chemical system. 
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Figure 6: Identification of a single time dependent source from instant measure- 
ments in one dimension with 1% noise in the data. Left: forward problem solu- 
tion u{x, t). Right: imaging functional J{s) for s = (x, t) G [0, 1] x [0, 0.2]. Hori- 
zontal axis is X, vertical axis is t. True source location is yellow x, measurement 
locations are yellow o, source position estimated by Algorithm [T] is black □. True 
source parameters (a,?/, r) are (3,0.4,0.05), estimated are (3.178,0.392,0.045). 



4.5.1 One dimensional case 



Let us consider a scalar forward problem of the form (|2.1|) with e 
L = 5, Q{u) 



1, ii; = 0, 



[0, 1] and T = 0.2. A single source of the form 
/(x, t) = aS{t — r)S{x — y) 



(4.4) 



is to be determined, where a = 3, r = 0.05 and y = 0.4. 

Similarly to the two dimensional case we use a finite difference scheme in 
space and an exponential integrator in time. To avoid committing an inverse 
crime we use a fine grid to compute the forward solution for simulating the data 
with 200 grid steps in x and 120 time steps. A coarser grid is used in Algorithm 
[l]to identify the source with 101 steps in both spatial and temporal variables. 

We observed from our numerical experiments that identification of time de- 
pendent sources is more sensitive to noise and numerical errors than the iden- 



tification of sources in examples in sections 4.3 and 4.4 Thus, for stable source 
identification we need to use more measurements than is required to just make 



the system (2.20) formally determined. The source in (4.4) is determined by 



three parameters, but we use six measurements for our numerical example. We 
make measurements at spatial locations 0.2, 0.5 and 0.8 at two time instants 
0.1 and 0.15 for a total of 6 measurements. 

The results of the forward simulation and source identification by Algorithm 
[l] are shown in Figure [6j We observe that both the source location and its 
intensity were identified reasonably well given the noisy data. The plot of the 
imaging functional J{s) explains why the time dependent source identification 



24 




0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 



0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 



Figure 7: Identification of a single time dependent source from instant measure- 
ments in two dimensions with 1% noise in the data. Left to right: shces of the 
imaging functional J(s), s = (y, t) for the three values of t = 0.008, 0.009, 0.010 
and y e ft. True source spatial location is yellow x, measurement spatial lo- 
cations are blue o. Source spatial position estimated by Algorithm [l] is black 
□ . True source parameters (a,x,?/,r) are (10,0.7,0.3,0.010), estimated are 
(11.18,0.73,0.31,0.009). 



problem is more difficult than the previously considered examples. The imag- 
ing functional has a large plateau surrounding the true source location, which 
decreases the discriminatory power of the method. For higher noise levels we 
observed that the estimated source location ends up somewhere on this plateau 
far away from the true source position. 

4.5.2 Two dimensional case 

Here we present the results of time dependent source identification for the three 
component chemical system in two dimensions from measurements instant in 
time. As we observed in Section 14.5.11 such source detection is more difficult 



than the cases considered in sections 4.3 and 14.41 so for stable identification 



we reduced the non-linearity of the system by taking smaller reaction rates 
ki = 100 and k2 = 200. Higher reaction rates leading to stiffer system can 
be handled using more efficient numerical schemes, for example [7 . However, 
proper numerical treatment of stiff systems is out of the scope of this work, so 
for convenience we work with reduced reaction rates in this section. 

We simulate the system up to T = 0.03, which is the time when the system 
is still in transient behavior. The source goes off at r = 0.01 and we make two 
sets of three point measurements each at instants = 0.015 and = 0.020 for 
a total of six measurements. 

In Figure [t] we show three slices of the imaging functional J{s) at time 
instants adjacent to the temporal source location estimated by Algorithm [l] We 
observe that the imaging functional has a narrow ridge with a plateau on top, 
which can make source identification difficult similarly to the one dimensional 
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case, where in the presence of noise Algorithm [T] can get stuck far away from 
the true source location. 

5 Conclusions and future work 

We presented here a method for source identification in non-linear time de- 
pendent advection-diffusion-reaction systems. We also provided the results of 
extensive numerical experiments that suggest that our method performs well in 
the presence of noise in the data and/or uncertainty in the number of sources 
present. The numerical experiments also show that the method's performance 
can be further improved by adaptively adding more measurements using the 
proposed strategy. 

The following topics of future study can be proposed. First, determining the 
conditions under which Algorithm [l] converges and proving the convergence. 
The analysis is complicated by the lack of regularity of solutions in the presence 
of point sources and by the coupling between the forward iteration and source 
estimation at each step of the algorithm. 

Second, the study of the case where only a partial knowledge of domain ft is 
assumed. In this case both the sources and the obstacles need to be determined. 
A method proposed in [3] solves the linear case by using the comparison results 
for elliptic equations. Since similar results hold for non-linear parabolic systems, 
it should be possible to extend the method in [3] to the setting considered here. 

Third, in this paper we assumed that all system parameters such as reac- 
tion rates, advection field and diffusion coefficients are known. In reality these 
parameters are estimated from some other measurements and thus are prone 
to inaccuracies. One may study the sensitivity of source identification with re- 
spect to uncertainties in the system parameters, or even try to estimate these 
parameters as a part of source identification problem. 
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